clear all; close all; clc
 addpath('../LibForSerial/LibrarySLP','../LibForSerial/LibraryLeSage/','../LibForSerial/LibVAR/');
%----------- Load shocks using Choleski identification ----------------%
  NumberConsecutiveshocks = 2; % 1 previous shocks that is  positive, and the 2nd one is positive too 
cd ./DataShocks/
  LoadShocksMandQpositiveEPUrevisedMAC; 
cd ../  
 sample=[ 1986 1; 2019  12]; ss_ = find(dates_EPU ==198601);
                             se_ = find(dates_EPU ==201912);
  uEPU   =   uEPU(ss_:se_);
  logEPU = logEPU(ss_:se_);
  IndPreviousNShocksEPU = IndPreviousNShocksEPU(ss_:se_);
  clear ss_ se_

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data = xlsread('../Dataset/VARdataBBupdatedMonthly.xlsx','MainData');
startsample=find(data(:,1)==sample(1,1)*100+sample(1,2));
endsample  =find(data(:,1)==sample(2,1)*100+sample(2,2));

data = data(startsample:endsample,2:end);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DATA TRANSFORMATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data(:,[3:10]) = log(data(:,[3:10]))*100;
%--------------------------------------%
clear startsample endsample 

T = size(data,2);

%P =   3;
 P  =  6; % including in the LP regression P lags of all variables in the system (3 lags in original BBD QJE; here 6 for consistency througout)
%P  = 12; % 
 
 plottype = 'Output';      ii=3; %
  
%% IRF 
h1 = 0 ;  %start LP at h=0 or h=1 (h=1 if impose no contemporaneous impact)
H  = 12*3; %#horizon of IR

% We are interested in the estimation of the dynamic multiplier of y_t+h with respect to a change in x_t for h ranging from 0 to H
%------------ IRF of GDP growth to a shock variable ----------------------%
    y  = data(:,ii);    %endogenous variable
%-------------------------------------------------------------------------%
  Remove_Large_Shocks = 0;
    x  = uEPU; 
    IndSeqShocks = IndPreviousNShocksEPU;
    IndLargeShocks = uEPU>2;

LargeShocksToBeRemoved = IndSeqShocks & IndLargeShocks;  
% [x IndSeqShocks IndLargeShocks]
if Remove_Large_Shocks; IndSeqShocks = IndSeqShocks - LargeShocksToBeRemoved; end
fprintf('Number of events with %d consecutive shocks is %d \n',NumberConsecutiveshocks,sum(IndSeqShocks));

%/////////////////////////////////////////////////////////////////////////%
%-------------------------------------------------------------------------%
xs = x.*IndSeqShocks; %shock variable times states
%-------------------------------------------------------------------------%
 
% w  = [lagmatrix( data(:,unique([ii   1 10])  ), 1:P )                       ]; % we control for employment as in BBD QJE VAR; otherwise IP response to 1 shock is positive (as per JMN AEJ)!
  w  = [lagmatrix( data(:,unique([ii 6 1 10 3])), 1:P ) lagmatrix(logEPU,1:P )]; %                     Replicating BBD QJE VAR system: SP500, FFR, EMP, IP, and EPU

%/////////////////////////////////////////////////////////////////////////%
%Additional robustness not reported in the paper; RUN with benchmark case P=6;
%/////////////////////////////////////////////////////////////////////////%
%%w  = [lagmatrix( data(:,unique([ii   1 10])  ), 1:P ) lagmatrix(logEPU,1:P )];
%/////////////////////////////////////////////////////////////////////////%
% The next 3 regressions study the effect of excluding employment as a control; IMPORTANT: the state multiplier is ALWAYS negative; 
% if anything             2 shocks confirm that uncertainty reduces activity, whereas 
% conclusions obtained by 1 shock are more sensitive to in/exclusion of employment
%%w  = [lagmatrix( data(:,unique([ii   1])     ), 1:P )                         ]; %control variables (contemporaneous vars, lagged vars)
%%w  = [lagmatrix( data(:,unique([ii 6 1  ])   ), 1:P )                         ];
%%w  = [lagmatrix( data(:,unique([ii   1  ])   ), 1:P )       data(:,6)         ];


if ii == 9
  w  = [lagmatrix( data(:,unique([ii 6 1 10 3])), 1:P )];
end
  w( ~isfinite(w) ) = 0;
OKforReg = isfinite(x) & isfinite(y) & isfinite(xs);
for cc=1:size(w,2)
    OKforReg = OKforReg & isfinite(w(:,cc)) ;
end  
  x = x(OKforReg,:); 
  y = y(OKforReg,:); 
  w = w(OKforReg,:);
  xs=xs(OKforReg,:);
  
% This shrinks the estimated dynamic multiplier to a polynomial of order r - 1
r = 3; %(r-1)=order of the limit polynomial (so r=3 implies the IR is shrunk towards a quadratic polynomial)
lambda = 20; %value of penalty
%-------------------------------------------------------------------------%
slp_sd  = locprojStateDep(y,x                  ,xs,w,h1,H,'smooth',r,lambda); %IR from Smooth Local Projection
 lp_sd  = locprojStateDep(y,x                  ,xs,w,h1,H,'reg');             %IR from        Local Projection

%% Figure 4, Panel (c) State-dependent IR to a EPU shock
figure(23)
%-------------------------------------------------------------------------%
     fname = sprintf('./DataShocks/EPUBBD_LP_linear.mat'); %generated from estimatebaselineLP_EPUv2.mat
load(fname)
if       ii==1;             slp_linear = slpShortRate;   slp_linearUp = slpShortRateUP;   slp_linearDown = slpShortRateDown; 
elseif   ii==3;             slp_linear = slpOutput;      slp_linearUp = slpOutputUP;      slp_linearDown = slpOutputDown;
elseif   ii==6 || ii==7;    slp_linear = slpStockMkt;    slp_linearUp = slpStockMktUP;    slp_linearDown = slpStockMktDown;
elseif   ii==5 || ii==9;    slp_linear = slpInflation;   slp_linearUp = slpInflationUP;   slp_linearDown = slpInflationDown;
end
%-------------------------------------------------------------------------%
%-------------------------------------------------------------------------%
subplot(1,2,1)
  plot( 0:H ,slp_linear, 'k', 0:H, slp_sd.IR+slp_sd.IRstate, 'b--', 0:H, slp_sd.IR,'k--', 'LineWidth', 1.5 ); 
  title(['IR_{slp} when there are ' num2str(NumberConsecutiveshocks) ' Consecutive Positive Shocks']); %legend(         's_t=+1 (Previous shock(s)=1)','s_t=0','Location','Best')
hold on; plot(0:H , zeros(H+1,1) , '-k' , 'LineWidth' , 2 ); grid; xlim([0 H]);                         legend('Linear','s_t=+1 (Previous shock(s)=1)','s_t=0','Location','Best')
if ii==2 
    ylim([-2.8 1.6]);set(gca,'YTick',[-2.8:0.5:0.6],'FontSize',12);ylabel('Percent','FontSize',12);
end
xlim([0      H]);set(gca,'XTick',[0 3:3:H],'FontSize',12);

subplot(1,2,2)
grpyat=[(0:1:H)', slp_sd.Interval_up_state; (H:-1:0)' slp_sd.Interval_down_state(end:-1:1)]; patch(grpyat(:,1), grpyat(:,2), [0.7 0.7 0.7],'edgecolor', [0.7 0.7 0.7]); hold on
plot(0  :H , zeros(H+1,1), 'k-'); hold on 
plot(    0:1:H,   slp_sd.IRstate, 'b--','LineWidth', 1.5); hold on
title('State multiplier of State-dependent IR to uncertainty shock')
grid;
xlim([0      H]);set(gca,'XTick',[0 3:3:H ],'FontSize',12);